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ABSTRACT 

We present new analytical estimates of the large scale bias of neutral hydrogen (HI) We use a 
simple, non-parametric model which monotonically relates the total mass of a halo M to t with its 
HI mass Mm at zero redshift; for earlier times we assume limiting models for the f2m evolution 
consistent with the data presently available, as well as two main scenarios for the evolution of our 
A/hi - M t ot relation. We find that both the linear and the first nonlinear bias terms exhibit a strong 
evolution with redshift, regardless of the specific limiting model assumed for the HI density over time. 
These analytical predictions are then shown to be consistent with measurements performed on the 
Millennium Simulation. Additionally, we show that this strong bias evolution does not sensibly affect 
the measurement of the H I power spectrum. 

Subject headings: cosmology: large-scale structure of universe, diffuse radiation 



1. INTRODUCTION 

The 2 1 cm emission line of neutral hydrogen has been 
the workhorse of astronomy for over 50 years. Its 
use for galactic astronomy cannot be overestimated, 
and it now promises to revo l utionize cosmology as 



and it now promises to revo l utionize cosmology as 
well fcf. iFurlanetto et all 120061: iBarkana fc Loebl 120071: 
Pritchard fc Loebl 120081: iMao et all 120081: IVisbal ct al. 



20081 and references therein) 



Compared to present probes of large-sc ale structure, 
the 21 cm intensity mappin g technique (jChang et al.1 
I2008at IWvithe et at] 12008a!) offers three advantages. 
First, it allows to reconstruct the spatial distribution of 
HI in the universe over volumes much larger than the 
ones currently probed by galaxy redshift surveys. Since 
the precision with which the power spectrum of density 
fluctuations can be measured is approximately propor- 
tional to the number of independent Fourier modes that 
fit into the survey volume, a larger survey volume leads 
to higher precision in the measur ed power spectrum an d 
hence cosmological parameters (jLoeb fc Wvithd I2008D . 
Second, 21 cm intensity mapping allows to probe a very 
wide range of rcdshifts and to map HI distribution in 
space and time not only for the epoch after reioniza- 
tion {z < 6) but also during reionization and even in 
the so-called dark ages, before the first galaxies were 
formed (z > 20). Third, using radio telescopes is one 
of the most economical ways to probe the dark energy 
evolution among all current and planned experiments, 
although several crucial advances in calibration and fore- 
ground removal must be made for the intensity mapping 
approach to reach its stated goals. 

As a dark energy probe, 21 cm intensity mapping 
is primarily useful for measuring the baryon acous- 
tic oscillations (BAO) in the matter power spectrum 
(IWvithe et al.ll2008bl : IChang et al.ll2008bl : IWvithe et all 
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I2008al lAnsari et all [200l iSeo et al.l[20Tol ). which pro- 
vides important information on constraining cosmologi- 
cal models and parameters. Measurements of the scale 
of the BAO in the matter power spectrum as a func- 
tion of redshift, when combined with the true physical 
scale, i.e., the sound horizon scale measured from the 
cosmic microwave background (CMB), probe the angu- 
lar diameter distances and Hubble parameters to vari- 
ous redshifts, and therefore dark energy properties (see , 
for instan ce. IHu fc White! 119961: lEisenstein et ail 1 1998 1: 
Eisensteinll20q2l: iBlake fc Glazebrooldl2003t iLinderl booa 
Hu fc Haimanl 120031 : ISeo fc Eisensteinl 120031 ). The dis- 
tinct oscillatory feature of the BAO allows us to iso- 
late and measure the BAO information independent of 
the overall broadband shape of the HI power spectrum. 
Therefore in principle BAO measurements do not require 
a knowledge of the HI bias, i.e., the relation between the 
neutral hydrogen and the total matter clustering, which 
affects the overall shape of the HI power spectrum. 

However, understanding the HI bias is important for 
the 21cm BAO surveys in the following ways. First the 
prediction of the signal to noise of the BAO measure- 
ments is sensitive to the knowledge of the c lustering bias 
of HI (|Chang et all l2008bt ISeo et all [2010h . Second, un- 
derstanding HI bias will help us to extract cosmologi- 
cal information from the broadband shape of the power 
spectrum in addition to the BAO. Third, the nonlin- 
ear effects on the BAO, such as the degradation and 
the shift of the feature, may differ for different biased 
tracers (e.g.. lPadmanabhan fc Whitei r2009: Me tha et al.1 
|2010| ) . Knowing the properties of HI bias will allow us to 
approximately estimate the degree of the expected non- 
linear effect on the BAO in the HI distribution. 

In this work, we focus on two critical aspects of the 
large-scale HI bias: its dependence on the relation be- 
tween H I mass and halo mass and its evolution with red- 
shift. To study the first aspect, we use a non-parametric 
model where we assume a one-to-one correspondence be- 
tween a H I mass present in a dark matter halo Mh i and 
the total mass M to t of this halo. Using this Mm — M to t 
relation we estimate the large-scale bias of the HI us- 
ing a model containing elements on the halo occupation 
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distribution (HOD) form alism and the HI mass func- 
tion measured at z — (IZwaan et al J 1200 5). We then 
use thi s relation to paint H I on Millennium Simulation 
halos (jSpringel et al.l 120051 ) and to measure the power 
spectrum of neutral hydrogen in order to test our model 
assumptions. 

With respect to the redshift evolution of the bias, it 
is absolutely critical to know how the H I mass function 
evolves with time. However, to the best of our knowl- 
edge, little is known about the redshift evolution of HI 
bo th on the theoretica l and the observational sides (but 
see lWvithe et a l. 2009). A detailed study of the bias evo- 
lution requires two extra pieces of information which are 
currently poorly constrained by observational data: the 
evolution of the total HI mass in the universe and the 
evolution of the HI mass function. For both of these, 
we empirically assume limiting evolution models consis- 
tent with the few available data that should bracket the 
actual evolution. We then proceed to predict the bias 
evolution in these scenarios, being aware of the fact that 
the actual evolution will lie somewhere in between these 
limiting cases. 

The paper is organized as follows. In we intro- 
duce the basics of our model, we derive the relation 
Mhi — Mtot at z = 0, we introduce the limiting cases 
for its evolution and we define the bias. In f]3]we present 
our main results: analytic predictions for the evolution 
of HI bias, measurement of its scale dependence on the 
Millennium Simulation and how it would affects mea- 
surements of the power spectrum. Finally, we conclude 
in 21 To compare analytical and numerical results, 
throughout the paper we use a flat ACDM cosmology 
with (Im = 0.25, as = 0.9, Hubble parameter h = 0.73 
and initial power spectrum index of n s = 1, consistent 
with the cosmological parameters used in the Millennium 
Simulation. 

2. BACKGROUND 

2.1. Basic formalism 

The HOD mod el (cf. ICoorav fc Shethl l200l 

IBerlind fc Weinberg] I2002L and references therein) 
provides a parameterized prescription for the spatial 
distribution of objects that can be found inside dark 
matter halos. Originally inspired to model the distribu- 
tion of galaxies, it can also be used to model the spatial 
distributions of other astronomical objects, such as the 
neutral hydrogen clouds. The advantage in using this 
approach is that the problem of deriving the spatial 
distribution naturally falls into two separate parts: the 
first part is concerned with the distribution of halos in 
the universe, and the second part involves how objects 
(i.e., galaxies, hydrogen clouds, etc) populate these 
halos. In this work we use some of the elements of this 
HOD formalism to attempt to estimate the large-scale 
bias, i.e. in a quasi-linear regime, where the estimations 
of the 21 cm correlations for studying the BAO will take 
place. 

The first component is a model for the distribution of 
dark halos in the universe. In the analytic approach, 
one specifies the halo mass function, n(M), the spatial 
correlations of the halos, and the density profiles for halos 
of mass M, all based on halo collapse models and fits to 
A-body simulations. 



To describe the mass function of dark matter halos (the 
spatial density of halos as a functio n of mass), we use 
the Sheth fc Tormen mass function (|Sheth et al.l 120011 : 
ISheth fc Tor7ner]l200l : 
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where p m = fl m Pc is the average matter density in the 
universe at the epoch of observation, v(M), defined by 
v(M) = S c /a(M,z), is the ratio between the threshold 
for growth of linear fluctuations 5 C = 1.69 and the rms 
variance of a sphere of radius R (corresponding to a mass 
M) at redshift z 
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Here, P(k, z) is the linear matter power spectrum and 
W(k 7 R) is the Fourier transform of a top-hat window of 
radius R. The f unction f{v) motivated by the ellipsoidal 
collapse model (jSheth fc Tormenll2002D is 



1 2av 2 r , , N 
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(3) 



where for normalization purposes, A = 0.322, p = 0.3 
and a = 0.707. 

The second component is a model specifying how the 
objects under consideration (in the present case, the 
amount of neutral hydrogen) occupy the dark halos as 
a function of halo mass. In a HOD model, for the largest 
scales, we would be interested in the average number 
of galaxies (A(M to t)) that populate a particular halo of 
mass M to t- In our case, we arc interested in the average 
number of HI atoms, or equivalently the average mass 
of HI inside a halo. For the scope of this paper, which 
deals with large scale-bias, we are not using the complete 
HOD prescription: this would mean specifying the com- 
plete probability distribution function of P(Mni\M tot ), 
where its higher moments contribute to the estimation of 
the bias on small, nonlinear scales; for large scale bias, we 
only need the average HI mass (Mjj i (Mtot)) ■ When in- 
terested in the small-scale bias, details the spatial and ve- 
locity distribution of the H I inside the halos are needed; 
these are out of the scope of this paper (see for in- 
stance iBagla fc Khandall 120091 lObreschkow et all 120091 
IWvithe et al.l 120091 for treatments that include small- 
scale clustering). 

In the following subsection, we explore different models 
for the neutral hydrogen mass function and the relations 
between the neutral hydrogen Mm and the total mass 
of a halo M to t. 

2.2. Neutral Hydrogen Mass Function 

The neutral hydrogen mass function of galaxies in the 
local universe has b een measured by the HIPASS survey 
(|Zwaan et al.1 12005). However, the masses of the halos 
that host galaxies with a given HI mass are not, gen- 
erally, known. Thus, we need to construct a plausible 
model for the relation between the total mass of a halo 
Mtot and its neutral hydrogen mass Mh i . 

While such a relation does not have to be a sim- 
ple function, models that assume a one-to-one cor- 
respondence between the galaxy luminosity and the 
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Fig. 1. — Fraction of blue galaxies as a function of total halo 
mass from matching the cumulative halo mass function and the 
observed SDSS galaxy luminosity function. The horizontal range 
shows the observed fraction on blue BCGs in the maxBCG galaxy 
cluster catalog (E. Rykoff, private communication). 



halo mass (IColm et all [19991 IKravtsov fe Klvpinl Il99l 
iKravtsov et al.l 12004ft provide remarkably good fit not 
only t o the galaxy correl ation functions (jConrov et al.1 
120061: iMarin et al] 120081) bu t also to a wide vari- 
ety of observ a tional tests dNagai fe Kravtsovl 120051: 
Gnedinl [20081 IVolonteri fe Gnedinl 120091: IGuo et al 



2009H G onrov fc Wechslerll2009HTinker fc Conrovl^oM T 
Therefore, it is tempting to assume such a relation be- 
tween the H I mass and total halo mass as well, that is, 
to assign one HIPASS source to one dark matter halo. 

There is one complication with such an approach - not 
all galaxies contain HI. Thus, we need to account for 
the fraction of galaxies that do. Unfortunately, there are 
no reliable observational measurements of the fraction of 
Hi-rich galaxies as a function of galaxy mass or lumi- 
nosity. Instead, we adopt a simple ansatz that all blue 
galaxies contain H I and all red galaxies do not. Then, we 
can use the abundance matching technique if we match 
the abundances of HIPASS and blue galaxies. 

Specifically, the relation between the total mass of a 
dark matter halo M to t and its H I mass Mh i can be ob- 
tained by matching the two cumulative mass functions, 

(4) 

where nnt is the observed cumulative mass function 
(jZwaan et alJl2005l ) and 



tihi(> M H i) = n B (> Mtot), 



n B (> M to t) 



dMf B (M)n h (M). 



(5) 



In the last equation, f B is the fraction of blue galaxies 
as a function of halo mass. 

That fraction can be computed from the mea- 
sured Sloan Digital Sky Surv ey (SDSS) red and blue 
galax y luminosity functions (jMontero-Dorta fc Pradal 
120091) using the same abundance-matching technique. 
Namely, using the total S PSS l uminosity function from 
iMontero-Dorta fc Pradal ([20091 ) we first match halo 
masses to r-band magnitudes as 



} ga i(< m r 



dMn h {M); 



(6) 



A/to 



such a matching gives us a relation between a galaxy 
magnitude m r and the halo mass M to t. Then, for each 
m r (and, hence, M to t) we compute the fraction of blue 
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Fig. 2. — Cumulative mass function of dark matter halos of blue 
galaxi es (top panel), the cumulative H I mass function of galaxies 
(from IZwaan et al.l 120051 . right panel) , and the relation between 
the galaxy HI mass and its halo mass (central panel) obtained by 
matching the two, at z = 0. The dotted line in the top panel shows 
the cumulative mass function of all dark matter halos, hosting both 
blue and red galaxies. 

galax ies from the measured blue gala xy luminosity func- 
tion ([Montero-Dorta fc Prad a 2009). The so- measured 
blue galaxy fraction is shown in Figure Q] together with 
am estimate of the blue BCG fraction from the maxBCG 
galaxy cluster catalog (E. Rykoff, private communica- 
tion). Estimates of the blue BCG fraction in X-ray se- 
lected clusters usua lly give substantially larger numbers 
(jBildfell et al.ll2008D : however, because X -ray emission is 
likely to correlate with the existence of a cooling flow 
(and, hence, a blue BCG), the latter estimates should 
be taken as upper limits (E. Rykoff, private communica- 
tion). 

The resulting matching between the HI mass of a 
galaxy and its total mass from Equation ([4]) is shown in 
Figure [2] The relation between Mh i and Mot from the 
central panel allows us to compute the clustering prop- 
erties of H I from the known clustering properties of the 
dark matter halos. 

Unfortunately, little is known about the H I mass func- 
tion beyond z ~ 0.1. To obtain the Mhi - M to t match- 
ing at higher redshifts, we therefore must extrapolate 
the observed HI mass function to earlier times. This ex- 
trapolation is, of course, not unique and it requires two 
ingredients: the evolution of the total HI mass and how 
the HI mass is distributed among halos. To make the 
matter even worse, the observational constraints on the 
total HI mass density in the universe phi (or, equiva- 
lently, the cosmological density parameter fini) are im- 
precise. Figure [3] shows a representat ive sample of these 
constraints for a range of redshifts (jZwaan et al.l 120051 : 
iProchaska et~^[2005TlRao et al.ll2006[ ). While at z « 
and z > 2 the measurements are reasonably robust, at 
intermediate redshifts the error bars are large and the 
behavior is non-monotonic. Therefore, to explore a wide 
enough range of possible evolutionary histories of Ohi ; 
we choose three representative models. 

• Model A provides a reasonable fit to all data 



4 



Marin et al. 



1.5 - 



Model A 
Model B 
Model C 




Fig. 3. — Three models (A, B, and C) for the redshift evolution of 
the total H I density that we explore in this paper. The lines for the 
models arc slightly shifted vertically for clarity - all models adopt 
f2jji = lO - ^ at z > 3. Black sy mbols with er r or ba rs show the 
observational m easurements from Zwaan ct_aL 



(2005) (triangle), 
(2005) (squares). 



I Rao c t al. (200f|) (circles), and Prochaska^^aT 



points except for a s ingle z 
iProchaska et all (|2005t) . 



2.3 point from 



• Model B follo ws the suggestion of 

IProchaska fc Wolfd (|2009fl . who argue that 
SIhi does not evolve between z = and z « 2. 
Since it provides a plausible lower limit to the 
redshift evolution of Sim, we us e this model eve n 
if it ignores the measurements of lRao et al.l (|2006[ ). 

• In the same spirit, Model C serves as an upper limit 
to the possible flm(z). Formally, it is inconsistent 
with the HIPASS measurement, but we assume 
that f^Hi in Model C increases rapidly from the 
HIPASS value of 3.8 x 1CT 4 to about 1(T 3 within 
a redshift interval Az C 1. 

To fully specify the HI mass function we then need 
to address how the total HI mass is distributed among 
halos as a function of redshift. If the shape of the HI 
mass function is approximately preserved at z > 0, then 

n H i(> A/hi, z) = q N (z)n {> 9 A /(z)Af H i), (7) 

where nn(> Mht) is the observed cumulative HIPASS 
mass function (| Zwaan et alJ l2005h at z = 0. The two 
quantities q^ and qM are not independent; for a given 
evolutionary history of neutral hydrogen contents of the 
universe pm{z) = £lm{z)p C rit, <1n and qM are related by 
the constraint 



Phi(^) : 



uMm 
nni(> Mm,z)dMm 



o 



PHl(0). 



(8) 



Obviously, q N (0) = Qm(0) = 1- 

It is not practical to explore all possible evolutionary 
histories for riHi(> Mm,z); therefore, we restrict this 
work to only two scenarios. 

Pure Number Density Evolution: scenario (here- 
after, PNE) assumes that the HI masses of 



individual halos do not change at high redshift 
(qm = 1 for all z), but only the number density 
of halos of a given mass changes to reflect the 
change in the total neutral hydrogen content of 
the universe, qw = ^Hi(z)/f2Hi(0). 

Pure Mass Evolution: scenario (hereafter, PME -a 
direct analog of a "pure luminosity evolution" sce- 
nario in optical surveys) adopts an opposite ex- 
treme where the number density of halos remain 
the same (q^ = 1 for all z), but, instead, their 
masses change, qm = ^hi(0)/Qhi(z)- 

The justification for such a restriction is as follows. If 
the number density of galaxies with a given H I mass in- 
creases monotically with increasing redshift, then qN{z) 
is a monotonic function of z (which is expected, since 
galaxies merge but do not break up). If, for a given halo, 
its HI also increases monotically with z (which is ex- 
pected, if HI gets consumed by star formation and new 
accretion never exceeds the consumption), then qM (z) is 
also a monotonic function of z. Since 9jv(0) = <7m(0) = 1, 
our two scenarios bracket all possible evolutionary path 
for gjv and qM- 

Obviously, these assumptions do not hold all the time 
-for example, for some galaxies accretion of fresh gas can 
exceed the consumption by star formation at some times. 
However, over long timescales and over the whole galaxy 
population, it appears plausible that our two scenarios 
bracket the majority of evolutionary histories for riHi(> 
M m ,z). 

One more ingredient of our model can evolve with red- 
shift -and that is the fraction of blue galaxies /s(M to t). 
The measured fraction of blue galaxies at z = from 
Figure [1] is exceptionally well (within the thickness of the 
line in Figure [T]) fitted by the following simple formula: 



f B (M tot ) = 1 + 



M t , 



Mi 



-0.25 



exp 



Aft, 



M, 



0.6N 



(9) 



with Mi = 2.5 x 1O 11 M and M 2 = 3.0 x 1O 14 M . It 
is tempting to associate Mi with the characteristic mass 
of dwarf ellipticals, and M 2 with the characteristic mass 
of a galaxy cluster. Then the behavior of Jb with M 
becomes well justified physically: for M tot > Mi a small 
fraction of all halos host dwarf and (at larger masses) 
regular elliptical galaxies, and fs becomes a gradually 
decreasing function of mass. But for M to t ^ M 2 halos 
turn into galaxy groups and clusters, and the fraction of 
halos dominated by blue galaxies plummets. 

With this plausible (although by no means proven) in- 
terpretation, we can come up with reasonable evolution- 
ary scenarios for fg as a function of z. Specifically, we 
adopt the parametric form (Equation ([9])), but allow Mi 
and Af 2 to evolve with redshift as 



M 1 (z) = 2.5 x 10 n (l4-z) a M Q , 
M 2 (z) = 3.0 x 10 14 (1 + z)^M Q . 



(10) 



It seems plausible to assume that the fraction of blue 
galaxies at a given halo mass is higher at high redshift, 
i.e., a > and j3 > 0. To investigate the dependence of 
our results on these two parameters, we consider below 
three representative cases with (a, /?) = (0, 0), (0, 1), and 
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(1, 1). This set of values is, of course, not exclusive, but 
it is sufficient to evaluate the robustness of our results to 
these parameters. 

In this work we therefore consider a range of possi- 
bilities -our models A, B, and C for H,m(z) and two 
evolutionary scenarios (PNE and PME) -as sampling a 
plausible range of the redshift dependence of the H I mass 
function. It seems reasonable to speculate that the actual 
evolution should lie somewhere in between these different 
scenarios. 

2.3. Large-Scale Bias of HI 

On large scales, and not considering redshift distor- 
tions, the HI bias parameters provide the relations be- 
tween the HI and d ark matter correlation functions 
(|Pan fc Szapudil[2005l) 

£,HI ~&l,Hl£dm, (11) 
C*/«&l,HlCdn, + &2,HI&! )HI (£d2dm + P*™.), (12) 

where £ and £ represent the 2-point and 3-point correla- 
tion functions, respectively. Similar relations are found 
for quantities calculated in Fourier space (power spec- 
trum and bispectrum). 

To calculate the large-scale bias of neutral hydrogen, 
we start by considering first the bias of dark matter halos 
with respect to the total matter clustering. In the Sheth- 
Tormen formalism, the bias parameters are prescriptions 
obtained from the spherical collapse model. For halos of 
ma ss M, the linear and (firs t) nonlinear terms are given 
by (|Scoccimarro et afll2001l) 



b'i(M): 
b h 2 (M): 



1 + ei 
8 

21 

with coefficients given by 



Ei, 

■ Ei ) + 62 + -E-2 ) 



3) 



(13) 
(14) 

(15) 
(16) 
(17) 
(18) 



where £ = 1,2 represent the linear and first nonlinear 
terms, where their redshift dependence is encompassed in 
v = 5 c /a(M, z), which depends on the redshift evolution 
of the linear power spectrum. 

Using the Mh i — M to t relation obtained in the previous 
section, the HI bias from all halos with a mass greater 
than M m j n is given by the following integral: 

1 f°° 

bi,m = / dMn h (M)b?(M) <M H i(M)> . (19) 

Phi Ji/ min 

where pn i is the density of neutral hydrogen: 



ei = 


au 2 - 1 


S c ' 


(2 = 


av 2 (av 2 — 


s 2 


Ei = 


2p 


5 C (1 + (av< 


E 2 




E x 


Oc 



Phi = / dMn h (M) (Mhi(M)) 

'M m ,„ 



(20) 



In what follows we present the analytic results for the 
large-scale bias and its evolution in redshift. 



3. RESULTS 
3.1. Analytic Results 

We calculate the bias parameter s as a function of red- 
shift for all models described in ^2.21 In Figure [4] we 
show the results for the models where the galaxy blue 
fraction does not evolve with redshift, i.e., in Equation 
(fT0|). a = /3 = 0. In the left panel we present the linear 
(6i) bias redshift evolution; in the right panel, the nonlin- 
ear (62) bias evolution, with M m ; n = 10 4 /i -1 Mq, which, 
given the fact that it is not expected that low mass halos 
(below 10 7 /i _1 Mo) contain significant amounts of HI, is 
set to the value mentioned above as a conservative limit. 
The black thick line shows the evolution of bias in the 
PME scenario; as for the PNE scenario, the red thin solid 
line shows the results for model A, the magenta dashed 
line shows results for model B, and the blue dotted line 
does it for model C. 

The first feature that arises from these analytical re- 
sults is the strong dependence of the bias terms on the 
redshift for all the models considered. Given that fini 
does not evolve significantly with redshift, the reason for 
this behavior lies in the fact that preserving the match- 
ing relations (Equation at higher redshift has the 
consequence that the most massive halos, which are less 
numerous at high redshift (i.e., are more biased with re- 
spect to the dark matter distribution), carry more HI 
mass than the halos with the same M to t at z = 0. There- 
fore, it is expected for both b\ and 62 to show a strong 
redshift dependence. 

For the PNE scenario, we can see differences in both 
bias parameters, though they are relatively small. Be- 
tween models A and B, the differences are increasing 
with redshift, with a peak around z ~ 2 -at the mo- 
ment where f2ni(^) of model B is the lowest compared 
to the rest of the models; since the fraction of neutral 
hydrogen in that model is low, in the PNE scenario this 
implies fewer halos, and therefore the bias is expected to 
be larger. In model C, the difference in the bias evolu- 
tion compared to model A appears only at low redshifts, 
since it starts with a higher J7hi and therefore having a 
smaller bias on redshifts close to z = 0. The differences 
for the nonlinear bias term 62 are qualitatively similar 
but much milder. 

As opposed to the PNE scenario, there are no differ- 
ences between the models biases in the PME scenario; 
this is not surprising since, by definition, in that sce- 
nario just the masses of halos change, but not their total 
number density. Hence, irrespective of the nm(z) evo- 
lution, the behavior of the bias here depends only on the 
cosmological evolution. When we compare both scenar- 
ios, the differences in the bias parameter b\, albeit small 
in general, grow steadily with redshift reaching 10% level 
between z = 1 and 4; note here that there are differences 
between the PME scenario and the model B in the PNE 
scenario only after z ~ 2. For &2 between z = and 1 
there are almost no differences between the models, but 
after that the differences grow until reaching almost a 
factor of 2 by z — 4. In both cases the PME scenario has 
a larger bias; this difference is due to the fact that more 
halos are needed in the PNE scenario to account for the 
change of neutral hydrogen in the universe, which leads 
to a lower halo mass for a given HI mass and hence a 
smaller bias than in the PME scenario. 
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Fig. 4. — Bias coefficients of neutral hydrogen clouds as a function of redshift. Left: linear bias bi vs. redshift in the pure mass evolution 
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magenta line) and model C (blue dotted line). Right: first nonlinear bias term 62 vs. redshift for the same models described in the left 
plot. All coefficients are calculated with Af m j n = 10 4 h~ 1 Mq. 



As for the differences between the different blue frac- 
tion evolution models, and the case where we simply as- 
sume all galaxies contain HI, we show in Figure IXT1 the 
quantity 



.FID' 



(21) 



where i = 1,2, and the fiducial (FID) model corresponds 
to the model A in the PNE scenario, with (a, (3) = (0, 0). 
The black dotted line represents the /b(M, z) = 1, i.e. 
all galaxies carry HI; the blue short-dashed line repre- 
sents the evolution model with (a, (3) = (0, 1) , and the 
long-dashed line shows the results for the blue fraction 
evolution model with (a,/3) = (1, 1), both in the model 
A of the PNE scenario. For other models and scenarios 
the results are qualitatively similar: between the differ- 
ent blue fraction evolution models both the linear and 
nonlinear bias do not differ significantly. But when we 
compare our FID model to that where all galaxies carry 
neutral hydrogen, i.e., /s = 1, then the differences are 
more significant, considering a realistic blue galaxy frac- 
tion lowers the H I bias, since the neutral hydrogen avoids 
accumulating in the bigger halos and concentrates more 
in mid-size and small mass halos. 

The weak dependence of the HI bias on the model 
adopted for the neutral hydrogen evolution means that 
the intensity mapping observations will be able to pro- 
vide only a limited amount of information on the global 
evolution of the neutral hydrogen abundance in the uni- 
verse. On the other hand, from a cosmologist's point 
of view, this is good news: a detailed knowledge of the 
HI evolution may not be needed for studying the dark 
energy evolution in the radio intensity mapping experi- 
ments. Note that even though the differences in 62 be- 
tween the PME and PNE models at high redshift are 
larger than those in bi, measurements on higher-order 
correlations have larger uncertainties which would make 
the 62 estimations less reliable. Nevertheless, if good 
signal-to-noise measurements can be done, this param- 
eter can play a relevant role in determining the right 
evolution model. 



3.2. Linear Bias calculated from simulations 

In order to have a consistency check of our analyti- 
cal calculations, we turn to estim ate the linear HI bia s 
from the Millennium Simulation (jSpringel et al.ll2005D . 
The N-body simulation was carried out for 2I60 3 par- 
ticles with mass m p — 8.6 x I0 8 /i _1 M in a box with 
Lhox = 500 h^ 1 Mpc on the side, with the same cosmo- 
logical parameters we use for our analytical estimates. 

Following Equation (fTTj) . we can estimate the linear 
bias in the simulation by 



C S .i(fc) = 



' -Ffai(fc) 

Pdm(fc) : 



(22) 



where Pdm(k) is the real-space dark matter power spec- 
trum and Phi(^) is the real-space HI power spectrum. 
In addition to validate our analytical method, this mea- 
surement should allow us to estimate the scales below 
which nonlinear effects become relevant. 

We use the Mh i — M to t relation obtained in Sj2?2] to 
assign the neutral hydrogen mass to dark matter halos 
and calculate the HI power spectrum Phi- Note that 
due to the mass resolution of the simulation, our dark 
matter halo catalog does not resolve well the halos below 
a certain mass threshold, and those halos may contribute 
substantially to the value of the HI bias. Therefore, con- 
sistency requires the bias measured in the Millennium 
Simulation to be compared to the value obtained ana- 
lytically through Equation (|I9|) assuming the same mass 
threshold M m j n . For this reason, we compare analytical 
estimates and measured values of the bias obtained by 
assuming M m i„ ~ 10 n ft. -1 MQ, which should correspond 
to halos with N p ~ 100 particles. Another aspect when 
we are "painting" H I in the simulation is its distribution 
inside the halos. In our case, we place all hydrogen at 
the center of the halo, i.e we are not assuming diffuse HI 
emission inside. While this might not model accurately 
the small-scale clustering of the H I, what matters in the 
large-scale clustering are the details of the Mm — M tot re- 
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shifts from the Millennium Simulation (symbols) for the model A in 
the PNE scenario, and our analytical predictions for the linear bias 
kmi (horizontal lines). The dashed lines show analytical predic- 
tions for the whole set of halos (M m i n = 10* li -1 M0), while solid 
lines show the analytical calculations for the set of halos resolved 
in the Millennium Simulation (M min = 1O 11 /i _1 M ). 

lation; the spatial profile of the neutral hydrogen within 
a halo will have a very limited effect on the measured bias 
and its scale dependence on large scales we present in this 
paper (jCoorav fc Shethll200l [Schulz fc White!l2006t) . 

Figure [6] shows, in symbols, the measured bias b^f 1 (fc) 
for a scries of redshifts from z — to z — 3 for the model 
A in the PNE scenario. The solid lines represent &hi,i 
from our analytical modeling, with a minimum mass cut 
Af m in = 10 11 /i _1 M©; the dashed lines are the results of 
6 H i,i with M min =10 4 h -1 M . 

The bias measured from the simulation using Equation 
(|22|) should agree with &hi,i obtained from Equation (|19p 
on large scales; on smaller scales, a contribution from the 
clust ering within a halo results in a scale dependence in 
bias (ISeliakl 120001) which we refer to as a nonlinear bias 
effect. From the figure, we indeed observe a strong red- 
shift evolution of large-scale bias in agreement with the 
behavior found in our analytical calculations; the agree- 
ment is of the order of 5%. The small differences are due 



mainly to the fact that our analytical halo mass function, 
though is a good approximation, indeed has small differ- 
ences to the one measured in simulations. The figure 
also shows that the linear bias assumptions break down 
for k > 0.15 /iMpc -1 , but the effective nonlinear scale 
(where the linear bias deviates strongly from a constant 
value) appears to evolve little with redshift, likely a con- 
sequence of our adopted Mhi — M to t relationship, given 
the fact that in simulations this scale varies greatly with 
redshift. We point out that if our assumption of the one- 
to-one correspondence between the HIPASS source and 
the halo breaks down for large halo masses, it will alter 
6jjf to some extent. 

As stated earlier, the difference between the two sets of 
horizontal lines in Figure [5] shows the contribution of ha- 
los with M tot < M min = lO^h^MQ. This contribution 
is substantial. It is important to note that in the case of 
the 21 cm intensity mapping, we suffer no lower detection 
limit in neutral hydrogen mass and therefore halo mass. 
In that case, the bias derived with M m j n = is the ap- 
propriate one. As a much higher mass resolution would 
be required to reduce M min much below 10 11 Ii~ 1 Mq, the 
Millennium Simulation -or another simulation with sim- 
ilar resolution- cannot be used to provide an accurate 
estimate for the H I bias for this kind of survey, but only 
as a useful tool to validate the semi-analytical estimates 
obtained in the previous section. 

3.3. Impact of bias evolution on the power spectrum 

Having shown how the HI bias evolves with redshift, 
we turn to assess to what extent this may or may not have 
an effect on the measurement of the power spectrum. At 
tree-level in perturbation theory, the density fluctuations 
in the neutral hydrogen <5h i and in the dark matter ddm = 
S are related through 



6 K 1 (S, X) = b(x)S(x, X ) = b(x)D+ ( X )S(x, 0) 



(23) 
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where D + (x) is the linear growth factor, and we ex- 
changed rcdshift z for comoving distance x- When ob- 
servations are performed, only the product E(x) can 
be measured, as the redshift evolution of the bias con- 
tributes a component that combines with the growth of 
structure. To understand to what extent this has an in- 
fluence on the measurement of the power spectrum, we 

start by Taylor expanding E(x) = Y^-^o^X" where 

all derivatives Eq = d n E/dx n are evaluated at the ob- 
server's position. Next, we Fourier transform 5m(x) to 
obtain 



d 3 xe ik - S S H i(x,x) 



d 3 xe- ikS 



E 



Z7"( n ) 

#o X 



00 in pi n ) 

■■Y,-i3-w> )> 



(24) 



n=0 



•nnn e -ik-2 Ncxt 

K11 



where we used the fact that x n 
we note that S(k, 0) is a Gaussian random variable with 



dispersion er^g ^ 



y / P(k, 0). We can then rewrite it as 
S(k, 0) = \J P(k, 0)Aj with Ag a Gaussian random vari- 
able with (Ag) = and (AgAt ) = (2ir) 3 6 a D (k - k'). Note 
that k for Ag is nothing but a label and that all the de- 
pendence of 8{k, 0) on the wavevector is really encoded 
in the prefactor. Defining the HI power spectrum by 

(SHi(k,x)S* R1 (k', x )} =Pni(k, X ) (2tt) 3 5%(k- k') we can 
then express Phi as 

Pai(k, X ) « [L(k)+v 2 H(k)] P(k,0), (25) 

where [i = cos(#) is the cosine of the angle the wavevector 
k makes with the line of sight and L(k) and H (k) are the 
first two terms that encode the x dependence of <5hi- We 
calculate the first few terms of L and H using the fact 

ksm(8) 



that fc|| = fccos((9) = 
They turn out to be 

L 



rp2 771 rp'f 



kfi and k± 
P' 



2kP' 



(26) 



Hf 



/ p/2 pll pi \ pl2 

where P' and P" denote first and second derivatives of P 
with respect to k, while E' and E' Q ' denote derivatives of 
E(x) with respect to the comoving distance evaluated at 
the observer position. In general, the comoving distance 
dependence generates in the power spectrum a direction- 
ality dependence that is reminiscent of redshift space dis- 
tortions. Despite that the nature of the effect considered 
here is completely different from the one of redshift space 
distortions (which is due to peculiar velocities), the two 
effects may be in practice difficult to separate from one 
another, unless we have a sufficiently small sample vari- 
ance on the clustering measurement. However, we point 
out that the n-th derivative with respect to the comov- 
ing distance scales as and therefore E' ~ 10~ 4 and 
Eq ~ 10~ 8 . As such, all but the first term appearing 
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Fig. 7. — First two correction coefficients L(k) and \H(k)\. 

in Equations ([26]) and (|2T[) are of order 10 -8 and only 
the Eq term appearing in L will give a non-negligible 
contribution. Also, it is possible to note that the above 
treatment is valid even in the case of a perfect unbiased 
dark matter tracer. 

In Figure [7] we show the functions L(k) and \H(k)\. 
Consistent with the above estimates, H turns out to be 
extremely small. On the other hand, L is dominated 
by the term Eq and therefore does not show any appre- 
ciable k dependence, except for extremely large scales 
k ~ 10 -3 , where it deviates from Eq by about 0.3%. 
These two facts together allow us to conclude that the 
measurement of the 21 cm power spectrum through in- 
tensity mapping should not be sensibly affected by the 
redshift dependence of the bias. 

4. SUMMARY AND CONCLUSIONS 

In this paper we have estimated analytically and in the 
Millennium Simulation the large-scale bias of neutral hy- 
drogen from a non-parametric model which establish a 
one-to-one relation between the total mass and the HI 
mass of galactic halos. Although only an approxima- 
tion, this relation is expected to hold for a large range of 
masses. 

The bias parameters have been calculated using a 
model which uses elements of the HOD formalism, in a 
way that can be considered complementary to the galaxy- 
HOD model of HI galaxies carried out bv lWvithe et al.l 
(120011 : we do not deal in our case with the individual 
galaxies inside dark matter halos but with the total neu- 
tral hydrogen inside them. In terms of large-scale clus- 
tering, these approaches are completely equivalent . 

Using this approach together with different models for 
the redshift evolution of the HI content in the universe 
and the relation between neutral hydrogen and total mat- 
ter contents of galactic halos, we analyze the evolution 
of bias with redshift. We find that both linear and non- 
linear bias increase signific antly with red s hift. T his con- 
trasts with the finding by iWvithe et al.l (|2008bD , based 
on an extr apolation to lower redsh ifts from a reionization 
model of IWvithe k Loebl pOOl . Since IWvithe et ahl 
(|2008bft do not account for clustering of galactic halos 
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that host all neutral hydrogen after reionization, they 
underestimate the HI bias at z < 4. On the other hand, 
our results are in agreement with the high value of the 
linear lar ge-scale bias at redshift z ~ 3 measured in the 
model by IBagla fc Khandai ( 200^) where they populate 
H I in high-resolution N— body simulations. 

We also find that, given the limitations of our model at 
high masses, the estimation of the bias does not change 
significantly when the very high mass halos are neglected. 

Our main conclusion is that, despite a wide range of 
possible models for the evolution of HI mass function 
that we consider, in all models the bias evolution is sim- 
ilar -the neutral hydrogen distribution is mildly antibi- 
ased at z = but becomes strongly biased (&hi ~ 2) by 
z ~ 4. This result is encouraging for the planned radio 
intensity mapping experiments, since large bias implies a 
stronger 21 cm signal. Nevertheless, this strong redshift 
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